%% process_pde_stokes_t.m
function pde = process_pde_stokes_t(index)
switch index
    case 1
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y,t) (x.^2.*y.^2+exp(-y)).*cos(2*pi*t);
        pde.u1dx = @(x,y,t) x.*y.^2.*cos(t.*pi.*2.0).*2.0;
        pde.u1dy = @(x,y,t) -cos(t.*pi.*2.0).*(exp(-y)-x.^2.*y.*2.0);
        pde.u2 = @(x,y,t) (-(2/3)*x.*y.^3+2-pi*sin(pi*x)).*cos(2*pi*t);
        pde.u2dx = @(x,y,t) -cos(t.*pi.*2.0).*(pi.^2.*cos(x.*pi)+y.^3.*(2.0./3.0));
        pde.u2dy = @(x,y,t) x.*y.^2.*cos(t.*pi.*2.0).*-2.0;
        pde.p = @(x,y,t) -(2-pi*sin(pi*x)).*cos(2*pi*y).*cos(2*pi*t);
        pde.pdx = @(x,y,t) pi.^2.*cos(t.*pi.*2.0).*cos(x.*pi).*cos(y.*pi.*2.0);
        pde.pdy = @(x,y,t) pi.*cos(t.*pi.*2.0).*sin(y.*pi.*2.0).*(pi.*sin(x.*pi)-2.0).*-2.0;
        pde.f1 = @(x,y,t) -nu.*(cos(t.*pi.*2.0).*(exp(-y)+x.^2.*2.0)+y.^2.*cos(t.*pi.*2.0).*2.0)-pi.*sin(t.*pi.*2.0).*(exp(-y)+x.^2.*y.^2).*2.0+pi.^2.*cos(t.*pi.*2.0).*cos(x.*pi).*cos(y.*pi.*2.0);
        pde.f2 = @(x,y,t) pi.*sin(t.*pi.*2.0).*(pi.*sin(x.*pi)+x.*y.^3.*(2.0./3.0)-2.0).*2.0+nu.*cos(t.*pi.*2.0).*(x.*y.*4.0-pi.^3.*sin(x.*pi))-pi.*cos(t.*pi.*2.0).*sin(y.*pi.*2.0).*(pi.*sin(x.*pi)-2.0).*2.0;
        pde.domain = [0, 1, -0.25, 0, 0, 1];
        pde.bdary = [0, 0, 0, 0];
    case 2
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y,t) pi*sin(2*pi*y).*sin(pi*x).^2.*sin(t);
        pde.u1dx = @(x,y,t) 2*pi^2*cos(pi*x).*sin(pi*x).*sin(2*pi*y).*sin(t);
        pde.u1dy = @(x,y,t) 2*pi^2*cos(2*pi*y).*sin(pi*x).^2.*sin(t);
        pde.u2 = @(x,y,t) -pi*sin(2*pi*x).*sin(pi*y).^2.*sin(t);
        pde.u2dx = @(x,y,t) -2*pi^2*cos(2*pi*x).*sin(pi*y).^2.*sin(t);
        pde.u2dy = @(x,y,t) -2*pi^2*cos(pi*y).*sin(2*pi*x).*sin(pi*y).*sin(t);
        pde.p = @(x,y,t) 10*cos(pi*x).*sin(pi*y).*sin(t);
        pde.pdx = @(x,y,t) -10*pi*sin(pi*x).*sin(pi*y).*sin(t);
        pde.pdy = @(x,y,t) 10*pi*cos(pi*x).*cos(pi*y).*sin(t);
        pde.f1 = @(x,y,t) pi*sin(pi*x).^2.*sin(2*pi*y).*cos(t)-nu*(2*pi^3*cos(pi*x).^2.*sin(2*pi*y).*sin(t)-6*pi^3*sin(pi*x).^2.*sin(2*pi*y).*sin(t))-10*pi*sin(pi*x).*sin(pi*y).*sin(t);
        pde.f2 = @(x,y,t) nu*(2*pi^3*cos(pi*y).^2.*sin(2*pi*x).*sin(t)-6*pi^3*sin(2*pi*x).*sin(pi*y).^2.*sin(t))-pi*sin(2*pi*x).*sin(pi*y).^2.*cos(t)+10*pi*cos(pi*x).*cos(pi*y).*sin(t);
        pde.domain = [0, 1, 0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
    case 3
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y,t) 10.*x.^2.*y.*cos(t).*(2.*y-1).*(x-1).^2.*(y-1);
        pde.u1dx = @(x,y,t) 20*x.*y.*cos(t).*(2*y-1).*(x-1).^2.*(y-1)+10*x.^2.*y.*cos(t).*(2*x-2).*(2*y-1).*(y-1);
        pde.u1dy = @(x,y,t) 10*x.^2.*cos(t).*(2*y-1).*(x-1).^2.*(y-1)+20*x.^2.*y.*cos(t).*(x-1).^2.*(y-1)+10*x.^2.*y.*cos(t).*(2*y-1).*(x-1).^2;
        pde.u2 = @(x,y,t) -10.*x.*y.^2.*cos(t).*(2.*x-1).*(x-1).*(y-1).^2;
        pde.u2dx = @(x,y,t) -10*y.^2.*cos(t).*(2*x-1).*(x-1).*(y-1).^2-20*x.*y.^2.*cos(t).*(x-1).*(y-1).^2-10*x.*y.^2.*cos(t).*(2*x-1).*(y-1).^2;
        pde.u2dy = @(x,y,t) -20.*x.*y.*cos(t).*(2*x-1).*(x-1).*(y-1).^2-10*x.*y.^2.*cos(t).*(2*x-1).*(2*y-2).*(x-1);
        pde.p = @(x,y,t) cos(t).*(20.*x-10).*(2.*y-1);
        pde.pdx = @(x,y,t) 20.*cos(t).*(2*y-1);
        pde.pdy = @(x,y,t) 2*cos(t).*(20*x-10);
        pde.f1 = @(x,y,t) 20*cos(t).*(2*y-1)-nu*(40*x.^2.*cos(t).*(x-1).^2.*(y-1)+20*x.^2.*cos(t).*(2*y-1).*(x-1).^2+40*x.^2.*y.*cos(t).*(x-1).^2+20*x.^2.*y.*cos(t).*(2*y-1).*(y-1)+20*y.*cos(t).*(2*y-1).*(x-1).^2.*(y-1)+40*x.*y.*cos(t).*(2*x-2).*(2*y-1).*(y-1))-10*x.^2.*y.*sin(t).*(2*y-1).*(x-1).^2.*(y-1);
        pde.f2 = @(x,y,t) 2*cos(t).*(20*x-10)+nu*(40*y.^2.*cos(t).*(x-1).*(y-1).^2+20*y.^2.*cos(t).*(2*x-1).*(y-1).^2+40*x.*y.^2.*cos(t).*(y-1).^2+20*x.*y.^2.*cos(t).*(2*x-1).*(x-1)+20*x.*cos(t).*(2*x-1).*(x-1).*(y-1).^2+40*x.*y.*cos(t).*(2*x-1).*(2*y-2).*(x-1))+10*x.*y.^2.*sin(t).*(2*x-1).*(x-1).*(y-1).^2;
        pde.domain = [0, 1, 0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
    case 4
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y,t) 5*(x.^2-2*x.^3+x.^4).*(2*y-6*y.^2+4*y.^3).*cos(2*pi*t);
        pde.u1dx = @(x,y,t) x.*y.*cos(t.*pi.*2.0).*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).*2.0e+1;
        pde.u1dy = @(x,y,t) x.^2.*cos(t.*pi.*2.0).*(x-1.0).^2.*(y.*-1.2e+1+y.^2.*1.2e+1+2.0).*5.0;
        pde.u2 = @(x,y,t) -5*(y.^2-2*y.^3+y.^4).*(2*x-6*x.^2+4*x.^3).*cos(2*pi*t);
        pde.u2dx = @(x,y,t) y.^2.*cos(t.*pi.*2.0).*(y-1.0).^2.*(x.*-1.2e+1+x.^2.*1.2e+1+2.0).*-5.0;
        pde.u2dy = @(x,y,t) x.*y.*cos(t.*pi.*2.0).*(x.*-3.0+x.^2.*2.0+1.0).*(y.*-3.0+y.^2.*2.0+1.0).*-2.0e+1;
        pde.p = @(x,y,t) 5*(2*x-1).*(2*y-1).*cos(2*pi*t);
        pde.pdx = @(x,y,t) cos(t.*pi.*2.0).*(y.*2.0-1.0).*1.0e+1;
        pde.pdy = @(x,y,t) cos(t.*pi.*2.0).*(x.*1.0e+1-5.0).*2.0;
        pde.f1 = @(x,y,t) cos(t.*pi.*2.0).*(y.*2.0-1.0).*1.0e+1-nu.*(cos(t.*pi.*2.0).*(y.*2.4e+1-1.2e+1).*(x.^2.*5.0-x.^3.*1.0e+1+x.^4.*5.0)+cos(t.*pi.*2.0).*(x.*-6.0e+1+x.^2.*6.0e+1+1.0e+1).*(y.*2.0-y.^2.*6.0+y.^3.*4.0))-pi.*sin(t.*pi.*2.0).*(x.^2.*5.0-x.^3.*1.0e+1+x.^4.*5.0).*(y.*2.0-y.^2.*6.0+y.^3.*4.0).*2.0;
        pde.f2 = @(x,y,t) cos(t.*pi.*2.0).*(x.*1.0e+1-5.0).*2.0+nu.*(cos(t.*pi.*2.0).*(x.*2.4e+1-1.2e+1).*(y.^2.*5.0-y.^3.*1.0e+1+y.^4.*5.0)+cos(t.*pi.*2.0).*(y.*-6.0e+1+y.^2.*6.0e+1+1.0e+1).*(x.*2.0-x.^2.*6.0+x.^3.*4.0))+pi.*sin(t.*pi.*2.0).*(y.^2.*5.0-y.^3.*1.0e+1+y.^4.*5.0).*(x.*2.0-x.^2.*6.0+x.^3.*4.0).*2.0;
        pde.domain = [0, 1, 0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
    case 5
        pde.index = index;
        nu = 1;
        pde.nu = nu;
        pde.u1 = @(x,y,t)-sin(pi*x).^2.*sin(2*pi*y).*cos(t);
        pde.u1dx = @(x,y,t)pi.*cos(x.*pi).*sin(x.*pi).*sin(y.*pi.*2.0).*cos(t).*-2.0;
        pde.u1dy = @(x,y,t)pi.*cos(y.*pi.*2.0).*sin(x.*pi).^2.*cos(t).*-2.0;
        pde.u2 = @(x,y,t)sin(2*pi*x).*sin(pi*y).^2.*cos(t);
        pde.u2dx = @(x,y,t)pi.*cos(x.*pi.*2.0).*sin(y.*pi).^2.*cos(t).*2.0;
        pde.u2dy = @(x,y,t)pi.*cos(y.*pi).*sin(x.*pi.*2.0).*sin(y.*pi).*cos(t).*2.0;
        pde.p = @(x,y,t)(x.*y-0.25).*cos(t);
        pde.pdx = @(x,y,t)y.*cos(t);
        pde.pdy = @(x,y,t)x.*cos(t);
        pde.f1 = @(x,y,t)nu.*(pi.^2.*cos(x.*pi).^2.*sin(y.*pi.*2.0).*cos(t).*2.0-pi.^2.*sin(x.*pi).^2.*sin(y.*pi.*2.0).*cos(t).*6.0)+y.*cos(t)+sin(x.*pi).^2.*sin(y.*pi.*2.0).*sin(t);
        pde.f2 = @(x,y,t)-nu.*(pi.^2.*cos(y.*pi).^2.*sin(x.*pi.*2.0).*cos(t).*2.0-pi.^2.*sin(x.*pi.*2.0).*sin(y.*pi).^2.*cos(t).*6.0)+x.*cos(t)-sin(x.*pi.*2.0).*sin(y.*pi).^2.*sin(t);
        pde.domain = [0, 1, 0, 1, 0, 1];
        pde.bdary = [0, 0, 0, 0];
    otherwise
        error("Invalid PDE index.");
end
fprintf("pde index: %d\n", pde.index);
fprintf("pde domain: xmin,xmax,ymin,ymax,tmin,tmax\n");
disp(pde.domain);
fprintf("pde boundary: bottom,right,top,left (0=Dirichlet,1=Neumann)\n");
disp(pde.bdary);
if isfield(pde, "bdary_p")
    fprintf("pde boundary of p: (0=Dirichlet,1=Neumann)\n");
    disp(pde.bdary_p);
end
end